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Abstract 

A  three-dimensional,  single-phase,  non-isothermal  numerical  model  for  proton  exchange  membrane  (PEM)  fuel  cell  at  high  operating  temperature 
(T  >  393  K)  was  developed  and  implemented  into  a  computational  fluid  dynamic  (CFD)  code.  The  model  accounts  for  convective  and  diffusive 
transport  and  allows  predicting  the  concentration  of  species.  The  heat  generated  from  electrochemical  reactions,  entropic  heat  and  ohmic  heat 
arising  from  the  electrolyte  ionic  resistance  were  considered.  The  heat  transport  model  was  coupled  with  the  electrochemical  and  mass  transport 
models.  The  product  water  was  assumed  to  be  vaporous  and  treated  as  ideal  gas.  Water  transportation  across  the  membrane  was  ignored  because 
of  its  low  water  electro-osmosis  drag  force  in  the  polymer  polybenzimidazole  (PBI)  membrane.  The  results  show  that  the  thermal  effects  strongly 
affect  the  fuel  cell  performance.  The  current  density  increases  with  the  increasing  of  operating  temperature.  In  addition,  numerical  prediction 
reveals  that  the  width  and  distribution  of  gas  channel  and  current  collector  land  area  are  key  optimization  parameters  for  the  cell  performance 
improvement. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Proton  exchange  membrane  (PEM)  fuel  cells  are  electro¬ 
chemical  devices  that  directly  convert  the  energy  from  the 
chemical  reaction  into  electricity.  Useful  features  such  as  high 
power  density,  simple,  safe  construction  and  fast  start-up  make 
those  particularly  suitable  for  home  appliance,  vehicles  and 
transportation  tools  [1].  Generally,  PEM  fuel  cells  operate  at 
temperature  below  363  K,  allowing  for  faster  start-up  and  imme¬ 
diate  response  to  changes  in  the  demand  for  power  [2] .  It  is  well 
known  that  the  operating  temperature  has  a  significant  influence 
on  PEM  fuel  cell  performance  [3].  The  increase  in  the  operat¬ 
ing  temperature  is  beneficial  to  fuel  cell  performance  since  it 
increases  reaction  rate  and  higher  mass  transfer  rate  but  usually 
lowers  cell  ohmic  resistance  arising  from  the  higher  ionic  con¬ 
ductivity  of  the  electrolyte  membrane.  In  addition,  at  high  tem¬ 
perature,  CO  poisoning  can  be  alleviated  by  reducing  chemisorp¬ 
tions  of  CO  [4].  A  great  deal  of  effort  has  been  expended  in  the 
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development  of  new  kinds  of  membranes  for  PEM  fuel  cells 
that  perform  efficiently  at  high  temperature  [5-10].  Recently, 
much  more  attention  has  been  paid  to  high  temperature  PEM 
fuel  cells-based  on  the  temperature  resistant  polymer  polybenz¬ 
imidazole  (PBI)  membrane,  which  allows  to  work  at  temperature 
up  to  473  K  [1 1,12].  Another  advantage  of  PBI-based  PEM  fuel 
cell  is  that  such  system  does  not  require  humidification  of  the 
membrane  [13];  traditional  membrane  such  as  sulfonated  flu- 
ourocarbon  polymer  (e.g.  Nafion)  requires  remaining  hydrated 
during  the  fuel  cell  reaction.  To  keep  the  membrane  hydrated, 
the  system  generally  includes  a  sub-system  that  humidifies  the 
cathode  air  stream  to  prevent  the  air  stream  from  drying  out 
the  membrane  as  the  air  is  flowed  through  the  fuel  cell.  There¬ 
fore,  development  of  PBI  membrane-based  fuel  cell  operating  at 
high  temperature  could  lead  a  simple  PEM  fuel  cell.  However, 
problematic  aspects  such  as  material  problems  related  to  cor¬ 
rosion,  electrode  degradation,  electrocatalyst  sintering  as  well 
as  re-crystallization  and  electrolyte  loss  by  evaporation  are  also 
accelerated  at  higher  operating  temperature.  These  material  con¬ 
straints  limit  the  temperatures  at  which  the  various  fuel  cells  can 
be  effectively  operated.  In  order  to  improve  and  optimize  the 
PEM  fuel  cell  design,  it  is  extremely  important  to  understand 
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the  thermal  and  electrochemical  behavior  under  various  design 
and  operating  conditions. 

The  experimentally  difficult  environment  of  fuel  cell  sys¬ 
tem  has  motivated  development  of  models  that  can  simulate  and 
predict  the  performance  of  PEM  fuel  cells.  Extensive  research 
efforts  have  been  devoted  in  the  past  decade.  The  first  one¬ 
dimensional  models  were  published  in  the  early  1990s  by 
Springer  et  al.  [14]  and  Bernardi  and  Verbrugge  [15].  More 
recently,  computational  fluid  dynamic  (CFD)  and  improved 
transport  models  have  made  the  development  of  more  realistic 
computational  model  feasible.  A  number  of  PEM  fuel  cell  mod¬ 
els  for  the  general  operating  temperature  cases  have  appeared 
in  the  literatures  [16-22].  However,  future  studies  are  still  in 
need,  especially  for  PEM  fuel  cell  operating  at  high  tempera¬ 
ture,  where  water  is  in  the  vapor  phase  and  the  membrane  has  a 
water  electro-osmosis  drag  force  near  zero.  Under  this  condition, 
the  performance  of  membrane  would  be  relatively  independent 
of  the  humidity,  and  thus  the  water  management  becomes  easier 
or  even  unnecessary  [23-25]. 

In  this  paper  a  single-phase,  three-dimensional,  non- 
isothermal  model  for  PEM  fuel  cell  at  high  operating  tem¬ 
perature  (T>  393  K)  is  presented  to  describe  the  fundamental 
processes  occurring  in  each  components  of  a  fuel  cell — current 
collector,  gas  flow  channels,  gas  diffusion  layers  (GDL),  cata¬ 
lyst  layers  (CL)  and  the  membrane.  Two  electric  potential  field 
equations  were  solved.  One  potential  field  was  solved  in  the 
membrane  and  catalyst  layers.  The  other  was  solved  in  the  cat¬ 
alyst  layers,  the  gas  diffusion  layers  and  the  current  collectors. 
The  Bulter-Volmer  equations  were  used  to  describe  the  rela¬ 


tionship  between  the  current  density  and  the  local  overpotential. 
The  convection  and  diffusion  of  different  species  in  the  gas  flow 
channel  and  gas  diffusion  layer  were  also  considered.  The  tem¬ 
perature  effects  on  mass  diffusivity  and  electric  conductivity 
were  taken  into  account.  The  heat  capacity,  gas  viscosity  and 
thermal  conductivities  of  each  gas  were  assumed  to  be  polyno¬ 
mial  functions  of  temperature.  The  model  was  implemented  into 
the  commercial  CFD  code  FLUENT  6. 1  with  custom  developed 
user-define  functions  (UDF)  [26]. 

2.  Physical  and  numerical  model 

Fig.  1  schematically  shows  a  PEM  fuel  cell  divided  into  the 
following  sub-regions:  the  current  collector,  gas  flow  channel, 
gas  diffusion  layer  and  catalyst  layer  in  the  anode  and  cathode 
sides  and  the  membrane  in  the  middle.  The  reactant  feed  is  con¬ 
veyed  by  the  gas  flow  channel  and  distributed  onto  the  anode  and 
cathode.  Reactants  pass  through  the  respective  porous  GDLs  and 
reach  the  CLs  where  the  electrochemical  reaction  occurs.  The 
membrane  acts  as  the  gas  separator,  electrolyte  and  the  proton 
conductor.  The  electrons  are  collected  by  the  anode  current  col¬ 
lector,  which  is  connected  to  cathode  current  collector  through 
the  external  load. 

2.7.  Assumptions 

In  this  model,  the  anode  feed  is  pure  hydrogen  and  air  is 
paralleled  in  the  cathode  gas  channel.  The  fuel  cell  is  assumed 
to  operate  in  steady  state  under  constant  load  conditions.  Since 
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Fig.  1.  Schematic  model  of  a  PEM  fuel  cell. 
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the  gas  streams  in  the  flow  channel  are  at  low  velocities  (or 
low  Reynolds  number),  laminar  flow  and  ideal  gas  behavior  are 
assumed.  Additional  assumptions  are  as  follows: 

(1)  All  water  produced  by  the  electrochemical  reaction  is 
assumed  to  be  vaporous  phase  due  to  the  high  operating 
temperature  and  water  transportation  across  the  membrane 
is  ignored  since  water  drag  coefficient  for  high  temperature 
membrane  is  low. 

(2)  Dilute  solution  theory  is  used  to  determine  the  species  dif¬ 
fusion. 

(3)  The  membrane  is  considered  impermeable  to  gases.  The 
crossover  of  reactant  gases  and  product  water  is  neglected. 

(4)  Ohmic  heating  in  the  current  collector,  GDL  and  CL  is 
neglected  because  of  their  high  electric  conductivity. 

(5)  Both  GDL  and  CL  are  assumed  to  be  homogeneous  and 
isotropic. 


empirical  expression  [27] 


TL75(l/Mk  +  1  /Mf!2 

m(Ei%)1/3+(EiV/3)2 


X  10“3, 


where  V\k  is  the  atomic  diffusion  volume,  and  the  value  of  ^  V\k 
is  given  by  Cussler  [27]. 

Energy  :  V  •  ( u(pE  +  p)) 


—  V  •  (^effvr  —  YhkJk  +  (reff  •  u)j  +  Sh,  (8) 

where  hk  is  the  enthalpy  of  species  k.  reff  is  the  effective  stress 
tensor,  which  can  be  ignored  due  to  the  low  velocity  of  laminar 

gas  flow.  Aeff  is  the  effective  thermal  conductivity  in  a  porous 
material  consisting  of  the  electrode  solid  matrix  and  gas,  which 
is  given  by 

7-eff  =  off  +  (1  —  s)XSi  (9) 


2.2.  Governing  equations 


A  steady  state,  single-phase,  isothermal  model  of  PEM  fuel 
cell  consists  of  five  principles  of  conservation:  mass,  momen¬ 
tum,  species,  energy  and  charge.  Thus,  the  governing  equations 
can  be  written,  in  vector  form,  as 


Continuity  :  V  •  (p«)  =  Sm, 


where  u  denotes  the  superficial  velocity  vector  in  the  porous 
media,  p  is  the  density  of  gas  mixture,  which  can  be  calculated 
by  using 


1 


where  pk  is  the  density  of  species  k  and  it  can  be  obtained  from 
the  ideal  gas  law  relation 


Pk 


pMk_ 

RT  ’ 


where  p  is  the  gas  pressure  and  M \  is  the  molecular  weight  and 
R  is  the  universal  gas  constant. 


1 

Momentum  :  V  •  ( puu )  =  —  S/p  +  V  •  r  +  Su,  (4) 

£Z 


where  s  is  the  porosity  of  the  electrode  materials,  r  is  the  viscous 
stress  tensor. 


Species:  V  •  (puYk)  =  V  •  Jk  +  Sk,  (5) 

where  Jk  is  the  diffusive  mass  flux  vector,  which  can  be  written 
as 


where  As  is  the  thermal  conductivity  of  the  electrode  solid  matrix 
and  Af  is  the  thermal  conductivity  of  the  gas,  which  can  be 
expressed  as  a  polynomial  function  of  temperature 

Af  =  A0  +  AiT  +  A2T2  +  A3T\  (10) 


where  A/,  /  =  0,  . . .,  3  can  be  determined  by  the  experiment  of 
real  gases,  as  seen  in  Appendix  A.  Similar  to  the  thermal  con¬ 
ductivity,  the  viscosity  and  heat  capacity  of  each  gas  species  can 
also  be  described  by  the  polynomial  expression  of  temperature 
with  empirical  coefficients. 

Charge  :  0  =  V  •  (/csoiV0soi)  +  Sso i,  (1 1) 

0  —  ^  ’  (k mcm  V 0mcm  )  T  ^mcm  ?  (12) 

where  solid  potential  Eq.  (11)  accounts  for  the  electron  transport 
through  the  electrode  solid  conductive  materials;  the  membrane 
potential  Eq.  (12)  represents  the  proton  transport  through  the 
membrane.  kso\  and  /cmem  are  electronic  conductivity  of  electrode 
and  ionic  conductivity  of  membrane. 

There  are  six  source  terms,  Sm,  Su,  Sk,  Sh ,  SSoi  and  SmQm, 
which  represent  various  volumetric  sources  or  sinks  arising  from 
each  sub-region  of  a  fuel  cell.  Details  of  the  various  source  terms 
are  summarized  in  Table  1.  It  shows  that  either  generation  or 
consumption  of  gas  species  k  and  creation  of  electric  current 
occurs  only  in  the  CL  where  the  electrochemical  reactions  take 
place.  Sm,  Sk,  Sh,  Sso\  and  SmQm  terms  are  therefore  related  to 
the  transfer  current  through  the  solid  conductive  materials  and 
the  membrane.  The  transfer  current  at  anode  and  cathode  can  be 
described  by  Butler- Volmer  equations  as  follows  [28]: 


•ref 

Ja  =  *a 


exp 


Uc  ph  \ 

RT  ) 


(13) 


N- 1 

Jk  =  -'ffpDkjVYj-  (6) 

j= i 

Here  Dkj  is  the  binary-diffusion  coefficient,  which  depends  on 
temperature  and  pressure  and  can  be  calculated  according  to  the 


jc 


X 


f 

V  RT 


(14) 
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Table  1 


Source  terms  for  continuity,  momentum,  species,  energy  and  charge  conservation  equations  in  various  regions  of  a  polymer  electrolyte  fuel  cell 


Source  terms 

Flow  channels 

GDL 

GCL 

Membrane 

Continuity 

o 

ii 

d 

5m  =  0 

5m  — 

■Sri  = 

Mho 

- ja ,  for  anode  side 

2  F 

Mo9  Mh9o 

- -  jc  H - jc,  for  cathode  side 

4  F  IF 

— 

Momentum 

o 

II 

Co 

s  _  ^ 

V  — 

ii 

^GDL 

ou  — 

kgcl 

Sk  = 

Mlh  •  V  TT 

"  ja,forH2 

2  F 

Mr>9  .  „  „ 

Species 

o 

II 

Co 

O 

II 

Co 

Mo  9 

Sk  =  — for02 
,,4F 
Mh9o 

sk  =  -^r  jo  torn 2° 
2F 


Energy 

o 

ii 

Co 

O 

II 

si 

CO 

5h  =  0,  for  anode  side 

1 2 

5h  =  - 

Jc 

5h  =  - T\As\  +  |  jc  ric  ,  for  cathode  side 

2  F 

%em 

5soi  =  0 

5soi  =  0 

5Soi  =  — Ja  on  anode  side 

5soi  =  0 

Change 

5mem  =  0 

5mem  =  0 

5soi  =  Jc  on  cathode  side 

5mem  =  Ja  on  anode  side 

5mem  =  — Jc  on  cathode  side 

5mem  =  0 

where  F  =  96,487  C  mol-1  is  the  Faraday  constant.  pu2 ,  po2  and 
/?h2o  are  the  partial  pressure  of  the  reactant  gases.  p°  is  the  stan¬ 
dard  pressure.  p^0  is  the  vapor  pressure  of  the  steam  at  the 
operating  temperature,  which  can  be  found  from  steam  Tables. 
tya  and  ac  are  the  anodic  and  cathodic  charge  transfer  coeffi¬ 
cients.  and  /£ef  is  the  reference  exchange  current  density, 
which  depends  on  the  local  temperature, 


water  is  in  gaseous  phase.  The  definition  of  au2 ,  ao2  and  <jh2o 
are 


P  h2  P  o2 

aU2  ~  pO  ’  a°2  ~~  pO  ’ 


au2o  = 


Pu2o 

Pu2o 


(20) 


The  temperature  dependence  of  the  membrane  conductivity 
can  be  accurately  described  by  the  Arrhenius  equation  [31] 


•ref 


,-ref 


a 


0exp 


(15) 


%em  —  K0  OXp 


Ea,k 

R 


(21) 


•ref 


,0  exP 


ea,c  ( i  n 

R  \T  To) 


(16) 


where  Ea,c  and  is  a,  a  is  the  active  energy  [29],  and  is 
the  anodic  and  cathodic  reference  exchange  current  densities 
at  reference  temperature  To,  see  in  Table  2.  p,  p \  and  p2  are 
empirically  determined  concentration  parameters  for  /3  =  0.25, 
Pi  =  0.5  and  p2  =  0.25  [28].  As  the  partial  pressure  decreases,  the 
exchange  current  density  also  decreases,  resulting  in  a  decrease 
of  cell  performance.  The  local  overpotential  r\ a  and  ijc  in  Eqs. 
(13)  and  (14)  can  be  expressed  as  [30] 


?7a  —  0sol  0mem>  (12) 

0c  =  0sol  0mem  Fq.  (18) 


Here  Vo  is  the  thermodynamic  equilibrium  potential,  which  can 
be  given  by  (T>  373.15  K) 

Vo  =  1-17  -  2.756  x  10 ~4(T  -  373.15) 


+  4.308  x  10“5  In 


a\\ 


(ao01/2 


«h2o 


(19) 


This  equilibrium  potential  is  calculated  from  thermodynamic 
data  of  reaction  enthalpy  and  entropy  changes  while  the  product 


where  Ea,k  is  the  activation  energy  and  k o  is  the  pre-exponential 
factor. 

The  heat  release  from  CL  of  the  PEM  fuel  cell  is  caused 
by  the  changes  of  enthalpy  and  irreversibility  related  to  charge 


Table  2 

Electrochemical  and  thermal  properties 


Parameter 

Symbol 

Value 

Unit 

Porosity  of  GDL 

e 

0.8 

— 

Porosity  of  GCL 

£ 

0.6 

— 

GDL/GCL  hydraulic  permeability 

1.0e-15 

— 

Membrane  ionic  conductivity 

kq 

12.99 

Sm-1 

GDL/GCL  electric  conductivity 

K 

103.3 

S  m-1 

Electrode  electric  conductivity 

K 

535 

Sm"1 

Anodic  charge  transfer  coefficient 

a  a 

1.0 

— 

Cathodic  charge  transfer  coefficient 

ac 

1.0 

— 

Anode  reference  exchange  current 

,-ref 

£a.O 

1.0e8 

Am-3 

density 

Cathode  reference  exchange  current 

,-ref 

C.O 

1.7e2 

Am  3 

density 

Thermal  conductivity  of  GCL 

1.7 

W  m-1  k-1 

Thermal  conductivity  of  membrane 

0.95 

W  m_1  k_1 

Thermal  conductivity  of  GDL 

1.7 

W  m_1  k_1 

Thermal  conductivity  of  current 

25 

W  m-1  k_1 

collector 

1186 


J.  Peng,  S.J.  Lee  /  Journal  of  Power  Sources  162  (2006)  1182-1191 


transfer  [21].  According  to  the  assumption,  the  ohmic  heating 
is  ignored  in  the  current  collector,  GDL  and  CL  but  considered 
in  the  membrane  due  to  the  relative  low  ionic  conductivity  of 
membrane.  Empirically,  the  change  of  entropy  As,  as  a  function 
of  temperature  T,  can  be  expressed  as  [2] 

As  =  33.64  +  4.52564  x  10 ~2T  -  2.98397  x  10~5T2 

+  3.40625  x  10"97’3  -  2.60417  x  lO~nT4.  (22) 

This  equation  is  a  good  approximate  for  temperature  from 
373  K  to  1 137  K. 

Once  membrane  potential,  </>mem  and  membrane  conductivity, 
/Cmem  are  obtained,  local  current  density,  /,  can  be  calculated  by 

J  —  VmemV0mem.  (23) 

The  average  current  density,  which  is  the  average  of  the  local 
current  density  over  the  entire  membrane,  can  be  obtained  by 

4vg  =  -t'—  f  I  ■  dA,  (24) 

Amem  J  Amem 

where  Amem  is  the  membrane  area. 

2.3.  Boundary  conditions 

On  the  inlet  boundaries  of  anode  and  cathode  flow  channels, 
the  stoichiometric  mass  flow  rate  and  mass  fractions  of  species 
were  prescribed  with  the  gas  temperature  equal  to  the  operating 
temperature,  as  seen  in  Table  3.  The  pressure  boundary  condi¬ 
tions  were  used  at  the  outlet.  As  there  is  no  protonic  current 
leaving  the  fuel  cell  through  any  external  boundary,  zero  flux 
boundary  condition  for  0mem  was  applied  on  all  outside  bound¬ 
aries.  However,  there  are  external  boundaries  on  the  anode  and 
cathode  side,  which  contact  with  the  external  electric  circuit  and 
the  electrical  current  is  generated  only  through  these  boundaries. 
Therefore,  a  prescribed  fixed  value  for  </>so i  was  used  on  these 
boundaries.  The  anode  side  was  set  to  zero  and  the  value  pre¬ 
scribed  on  the  cathode  side  is  the  cell  operating  voltage.  Zero 

Table  3 


Geometrical  and  operating  parameters 


Parameter 

Value 

Unit 

Cell  width 

3.4 

mm 

Channel  length 

235 

mm 

Channel  height 

0.7 

mm 

Anode  channel  width 

0.7 

mm 

Cathode  channel  width 

1.0 

mm 

Anode  GDL  thickness 

0.34 

mm 

Anode  GCL  thickness 

0.04 

mm 

Membrane  thickness 

0.065 

mm 

Cathode  GCL  thickness 

0.11 

mm 

Cathode  GDL  thickness 

0.34 

mm 

Electrode  height 

2.0 

mm 

Operating  temperature 

398-433 

K 

Anode  stoichiometric  mass  flow  rate 

1.5 

— 

Cathode  stoichiometric  mass  flow  rate 

2.0 

— 

Anode  outlet  pressure 

1.1 

atm 

Cathode  outlet  pressure 

1.1 

atm 

Anode  inlet  mass  fraction  H2 

1.0 

— 

Cathode  inlet  mass  fraction  02:N2 

0.22:0.78 

— 

flux  boundary  condition  for  0soi  was  applied  on  remaining  exter¬ 
nal  boundaries.  As  the  reactant  gases  are  dielectric,  zero  flux  of 
0Soi  and  0mem  are  enforced  at  the  interface  between  the  electrode 
and  the  flow  channel. 

Due  to  the  structure  of  the  FLUENT  CFD  code,  the  inter¬ 
face  between  the  membrane  and  CL  is  defined  as  a  wall  in 
order  to  prevent  any  crossover  of  species  through  the  mem¬ 
brane.  The  wall  has  a  fluid  region  on  each  side.  This  was 
implemented  by  creating  a  “shadow”  of  the  wall  cell  layer  by 
the  FLUENT  automatically.  According  to  the  physical  property 
of  electrolyte  membrane,  only  protons  are  allowed  to  transfer 
though  the  electrolyte  membrane.  Therefore,  zero  flux  of  0soi  is 
defined  at  this  internal  wall.  However,  for  the  0mem,  the  contin¬ 
uous  flux  boundary  conditions  should  be  applied  and  it  can  be 
expressed  as 

(k mem  ^ 0mem )  '  LI  |—  —  (jC- mem  ^ 0mem )  ’  LI  T  |  +  ,  (25) 

where  symbol  ‘  —  ’  and  ‘+’  stand  for  two  sides  of  the  internal 
wall,  n  denotes  the  exterior  normal  vector  of  the  internal  wall. 
The  flux  continuous  boundary  conditions  were  enforced  by  the 
custom  developed  UDF. 

On  the  interface  between  GDL  and  catalyst  layer,  zero  flux 
of  0mem  is  defined  because  protons  transfer  mainly  through 
the  electrolyte  materials.  However,  since  the  current  collector 
is  electric,  the  flux  of  (pso\  across  the  areas  should  be  continuous. 
The  UDF  is  implemented  to  enforce  the  continuous  flux  of  0soi 
and  it  can  be  expressed  as 

(^sol^0sol)  ’  LI  |  —  —  (/CsolV0sol)  ’  Ll~^~\-\~.  (26) 

As  shown  in  Fig.  1,  the  fuel  cell  with  single  straight  gas  chan¬ 
nel  was  employed,  and  the  fuel  cell  flow  plate  was  formed  with 
a  plurality  of  gas  channels.  Therefore,  the  symmetrical  bound¬ 
ary  condition  is  applied  on  the  side  walls  normal  to  v-axis  and 
no-slip  boundary  condition  is  applied  on  remaining  walls  for  the 
momentum  equations. 

The  model  was  implemented  via  a  set  of  user-defined  func¬ 
tions  in  a  commercial  CFD  code,  FLUENT  6.1,  which  is  a 
parallel  code  using  the  finite  volume  method  and  iterative  segre¬ 
gated  implicit  solver.  Second  order  discretization  schemes  were 


Fig.  2.  Polarization  curve:  comparison  of  simulations  and  experiments. 
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used  for  all  transport  equations.  More  details  of  the  numerical 
procedure  can  be  found  from  the  literature  of  former  researchers 
[32,33].  Stringent  numerical  tests  were  performed  to  ensure  that 
the  solutions  are  independent  of  the  grid  size.  A  mesh  with  about 
1 80,000  grid  points  was  found  to  provide  sufficient  spatial  res¬ 
olution.  The  solution  was  considered  to  be  convergent  when  the 
relative  error  in  each  field  between  two  consecutive  iterations 
was  less  than  10~6. 


3.  Results  and  discussion 

Fig.  2  shows  the  comparison  of  polarization  curves  obtained 
from  numerical  prediction  and  experimental  results  at  differ¬ 
ent  operating  temperatures.  In  both  cases,  the  cell  was  operated 
at  constant  temperature  and  without  any  humidification,  i.e. 
with  dry  gases.  It  can  be  seen  that  the  model  predicts  results 
close  to  the  experimental  data.  In  Eq.  (19),  it  can  be  seen  that 
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Fig.  3.  Distributions  for  current  density,  concentration  of  oxygen  and  local  overpotential  at  Vceii  =  0-6V,  /avg  =  0.485  A  cm-2  and  F=433K:  (a)  average  current 
density  distribution  at  the  membrane;  (b)  oxygen  molar  concentration  distribution  at  the  cathode  CL;  (c)  local  overpotential  r/a  distribution  at  anode  CL;  (d)  local 
overpotential  r/c  distribution  at  cathode  CL. 
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the  thermodynamic  equilibrium  potential  decreases  as  operat¬ 
ing  temperature  increases.  However,  with  the  increasing  current 
density,  the  cell  voltage  is  mainly  dominated  by  the  ohmic  loss. 
Membrane  ionic  conductivity  increases  with  the  increasing  of 
operating  temperature  (Eq.  (21)).  Therefore,  at  the  high  operat¬ 
ing  temperature,  ohmic  loss  decreases  and  the  cell  performance 
is  improved. 

3.1.  Isothermal  model 

Fig.  3  shows  the  average  current  density,  concentra¬ 
tion  and  local  overpotential  distribution  for  Vceii  =  0.6V, 
4vg  =  0.485  A  cm-2  obtained  using  the  isothermal  model.  Since 
the  electric  current  path  from  the  areas  of  CL  under  the  gas  chan¬ 
nel  is  longer  than  the  path  from  the  area  of  CL  under  the  current 
collector  land  areas,  the  current  density  maxima  locate  under  the 
current  collector  land  area  because  of  the  influence  of  the  ohmic 
losses  in  the  CL  and  GDL  (in  Fig.  3(a)).  The  concentration  of 
oxygen  under  the  current  collector  is  smaller  than  concentration 
under  the  gas  channel  due  to  consumption  (in  Fig.  3(b)).  The 
concentration  of  oxygen  in  the  cathode  CL  decreases  monoton- 
ically  along  with  the  gas  flow  direction  as  the  electrochemical 
reaction  proceeds.  Therefore,  the  current  density  also  decreases 
along  with  the  flow  direction. 

Fig.  3(c)  and  (d)  show  the  local  overpotential  (rja  and  0c)  dis¬ 
tribution  on  the  CL.  ija  decreases  along  with  the  flow  direction 
and  maxima  locate  under  the  current  collector  land  area.  This 
is  coincident  with  the  distribution  of  current  density.  However, 
for  the  local  overpotential  r]c  at  the  cathode  CL,  the  absolute 
value  increases  along  with  the  flow  direction  (seen  in  Fig.  3(d)). 
This  results  from  the  ohmic  loss  along  with  the  flow  direction 
decreased  by  the  decrease  in  current  density.  The  results  are  dif¬ 
ferent  from  those  in  [21],  where  only  the  overpotential  at  the 
cathode  is  considered  and  the  anodic  overpotential  is  assumed 
to  be  constant.  Fig.  4  shows  the  distribution  of  solid  potential 
0Soi  at  the  solid  electrode.  Because  of  the  lateral  electronic  resis¬ 
tance,  the  minima  solid  potential  locates  at  the  anode  CL  under 
the  gas  flow  channel,  while  the  maxima  appears  at  the  cath¬ 
ode  CL  under  the  gas  flow  channel,  resulting  in  the  relative 
slower  electrochemical  reactions  and  hence  lower  current  den¬ 
sity.  The  Aopotential  lines  are  normal  to  the  flow  channel  and 
side  walls  since  the  fuel  gas  is  assumed  to  be  insulated  and 
the  symmetrical  boundary  conditions  are  applied  on  the  side 
walls.  The  distribution  exhibits  gradient  in  both  v  and  y  direc¬ 
tion  due  to  the  non-uniform  local  current  production  in  the  CL 
and  shows  that  ohmic  losses  are  larger  in  the  CL  under  gas  chan¬ 
nels.  Fig.  5  shows  the  distribution  of  membrane  potential  0me m 
at  the  membrane.  As  indicated,  due  to  non-uniform  local  cur¬ 
rent  production  in  the  adjacent  CL,  the  gradients  in  both  v  and 
y  direction  also  exist.  The  membrane  potential  along  with  the 
anode  side  interface  is  not  uniform,  which  is  absolutely  different 
from  that  assumption  in  [21]. 

The  effect  of  the  width  of  gas  channel  on  cell  performance  is 
shown  in  Fig.  6  where  the  width  of  both  anode  and  cathode  gas 
channels  are  set  to  be  0.7  mm  and  1.0  mm.  The  average  current 
densities  are  /avg,o.7  =  0.526  A  cm-2  and  /aVg,i.o  =  0.474  A  cm-2 
with  operating  voltage  Vceii  =  0.6  V.  The  average  current  density 
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Fig.  4.  Solid  potential  distribution  in  the  current  collector,  GDL  and  CL  on  the 
(a)  anode  and  (b)  cathode  side  at  Vceii  =  0-6  V,  T=  433  K. 
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Fig.  5.  Membrane  potential  distribution  in  fuel  cell  membrane  at  Vceii  =  0-6  V, 
F=433  K. 


J.  Peng,  S.J.  Lee  /  Journal  of  Power  Sources  162  (2006)  1182-1191 


1189 


Fig.  6.  Current  density  distribution  in  the  lateral  direction  in  the  middle  of  the 
cell  with  different  channel  width  at  Vceii  =  0.6  V,  T  =  433  K. 

is  increased  with  the  decreasing  of  gas  channel.  The  current 
collector  land  area  is  enlarged  while  the  width  of  gas  chan¬ 
nel  is  decreased.  The  ohmic  losses  are  decreased  due  to  the 
fact  that  electric  current  path  from  CL  to  the  current  collec¬ 
tor  land  area  is  shortened.  Therefore,  the  average  current  den¬ 
sity  is  increased  at  the  operating  point  where  the  cell  voltage 
is  essentially  dominated  by  ohmic  loss.  On  the  other  hand, 
it  is  predictable  that  with  decreasing  of  gas  channel  width, 
the  average  current  density  will  be  decreased  in  the  operat¬ 
ing  point  where  the  cell  voltage  is  essentially  dominated  by 
mass  transport  of  fuel  gas.  From  Fig.  6,  it  can  be  seen  that 
the  maxima  values  of  current  density  locate  under  the  current 
collector  land.  Close  to  the  margin  region,  though  current  col¬ 


lector  land  area  with  anode  gas  channel  width  1 .0  mm  is  smaller 
than  that  with  anode  gas  channel  width  0.7  mm,  the  former 
current  density  is  larger  than  the  latter,  because  more  hydro¬ 
gen  gas  can  be  supplied  with  increasing  of  anode  gas  channel 
width. 

It  is  well  known  that  the  fuel  cell  itself  has  many  trade-off 
options.  The  pressure  loss  along  the  gas  channel  is  increased 
with  the  increasing  of  gas  channel  width  and  the  operating  para¬ 
sitic  power  consumption  is  increased.  Therefore,  it  is  critical  to 
optimize  the  width  and  distribution  of  gas  channel  and  current 
collector  land  area  in  order  to  improve  cell  performance  after 
the  cell  operating  point  is  determined. 

3.2.  N on-isothermal  model 

As  we  have  presented,  the  heat  is  released  from  CL  of  the 
PEM  fuel  cell  through  the  electrochemical  reaction  and  from 
the  membrane  by  ohmic  resistance.  Fig.  7  shows  the  tempera¬ 
ture  distribution  of  the  membrane  with  the  inlet  gas  temperature 
T=433K.  For  low  current  density  state,  (Fig.  7(a)),  the  tem¬ 
perature  difference  from  the  inlet  gas  temperature  is  small  and 
the  temperature  maxima  locate  under  the  gas  channel.  This  is 
caused  by  the  difference  of  thermal  conductivity  between  reac¬ 
tant  gases  and  the  electrode.  As  the  temperature  is  increased  with 
the  increase  in  current  density,  the  temperature  maxima  appear 
near  the  outlet  boundary  (Fig.  7(b)).  Fig.  8  shows  the  temper¬ 
ature  distribution  on  the  midway  section  of  the  fuel  cell.  The 
/^o-temperature  lines  are  normal  to  the  side  boundary  because 
symmetrical  boundary  conditions  were  applied.  In  low  current 
density  (Fig.  8(a)),  the  maxima  of  temperature  locate  in  the 
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Fig.  7.  Temperature  distribution  contour  on  the  membrane:  (a)  low  current  density  /avg  =  0.075  A  cm  2,  Vcen  =  0.8  V;  (b)  high  current  density  7av g  =  1.025  A  cm 
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Fig.  8.  Temperature  distribution  on  midway  section  of  the  fuel  cell:  (a) 
low  current  density  /avg  =  0.075  A  cm~2,  Vceii  =  0.8  V;  (b)  high  current  density 
/avg  =  1  -025  A  cm~2,  Vceu  =  0.4  V. 


cathode  CL  under  the  gas  channel.  However,  in  high  current 
density  (Fig.  8(b)),  it  is  observed  that  the  maxima  of  tempera¬ 
ture  shift  to  the  position  under  the  current  collector  land  area. 
This  can  be  explained  by  the  increased  ohmic  heat  in  the  mem¬ 
brane  at  high  current  density  state.  The  temperature  variation 
is  also  increased  with  the  increasing  of  current  density.  The 
influence  of  thermal  results  on  the  performance  of  fuel  cell  is 
shown  in  Fig.  9,  where  the  profiles  of  current  density  across 
the  membrane  are  presented.  From  Fig.  9(a),  the  current  den¬ 
sity  is  increased  when  the  heat  exchange  is  considered,  since  the 
temperature  inside  the  fuel  cell  is  larger  than  the  temperature 
at  the  surface  that  can  always  be  recognized  as  the  operating 
temperature.  The  difference  between  these  two  results  can  be 
seen  in  Fig.  9(b),  where  the  maxima  of  difference  appear  under 
the  collector  land  areas.  Therefore,  both  the  current  density  dis¬ 
tribution  and  temperature  distribution  are  closely  related  with 
the  geometry  structure  and  dimension  of  the  current  collector 
land  area. 
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Fig.  9.  The  effect  of  temperature  on  current  density  distribution  while 
Vceii  =  0-6  V:  (a)  current  density  distribution  with  isothermal  and  non-isothermal 
state;  (b)  relative  difference  between  the  computed  current  density  profiles  of 
isothermal  and  non-isothermal  state. 


4.  Conclusion 

In  this  paper,  a  three-dimensional,  single-phase,  non- 
isothermal  PEM  fuel  cell  model  at  high  operating  temperature 
(T>  393  K)  was  developed  and  implemented  in  the  framework 
of  a  CFD  code.  Water  was  considered  to  be  in  vaporous  phase 
and  the  water  transportation  across  the  membrane  was  ignored 
because  of  its  low  water  electro-osmosis  drag  force.  The  com¬ 
plete  set  of  conservation  equations  of  mass,  momentum,  energy, 
species  and  charge  were  numerically  solved  with  proper  account 
of  electrochemical  kinetics.  The  electron  transport  equation 
was  solved  in  the  CL,  GDL  and  current  collectors  rather  than 
assumed  uniform  and  constant,  rendering  more  accurate  predic¬ 
tion  of  local  overpotential  and  current  density. 

A  single  straight-channel  PEM  fuel  cell  at  operating  temper¬ 
ature  T=433K  was  numerically  studied  in  detail  with  focus  on 
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Table  A.  1 


Coefficients  of  polynomial  expression  of  gas  properties 


Item 

A() 

Ai 

a2 

^3 

Thermal  conductivity 

h2 

4.3484  x  1CT2 

4.8712  x  10“4 

-1.4917  x  10“7 

4.6636  x  10“ 11 

h2o 

5.1348  x  10“3 

1.8280  x  10“6 

1.5895  x  10“7 

-7.8887  x  10“n 

n2 

6.9408  x  10“4 

9.7432  x  10“5 

-5.0467  x  10"8 

1.5163  x  10“H 

o2 

5.5780  x  10“4 

9.5339  x  10“5 

-3.1382  x  10“8 

7.1220  x  10“12 

co2 

-8.5816  x  10"3 

8.5843  x  10“5 

-1.4887  x  10-9 

-7.9582  x  10“12 

Viscosity 

h2 

2.6488  x  10"6 

2.2381  x  10“8 

-5.0735  x  10"12 

8.2349  x  10“16 

h2o 

-3.1387  x  10-6 

4.1514  x  10“8 

0 

0 

n2 

3.3349  x  10“6 

5.4210  x  KT8 

-2.1159  x  10"11 

4.1614  x  10“15 

o2 

2.8879  x  10“6 

6.6299  x  10“8 

-2.5463  x  10“H 

4.9740  x  10“15 

co2 

7.5190  xl0“8 

5.5156  x  10“8 

-1.8831  x  10“n 

3.4136  x  10“15 

Heat  capacity 

h2 

1.3550  x  10-4 

3.6304 

-4.6474  x  10"3 

2.2471  x  10“6 

h2o 

2.0963  x  10“3 

-9.1474  x  10“‘ 

1.6648  x  10“3 

-5.5865  x  10“7 

n2 

1.0732  x  10-3 

-2.6917  x  10“' 

5.9522  x  10“4 

-2.3165  x  10"7 

o2 

8.3669  x  10-2 

2.5830  x  10“‘ 

7.7057  x  10“5 

-8.0416  x  10“8 

co2 

5.0437  x  10-2 

1.4174 

-9.0900  x  10"4 

2.2287  x  10“7 

the  temperature  distribution  and  cell  performance.  Overall,  the 
model  was  in  a  good  agreement  with  experimental  results.  The 
current  density  increases  with  the  increasing  of  operating  tem¬ 
perature,  demonstrating  a  necessity  for  non-isothermal  modeling 
of  PEM  fuel  cells.  The  maxima  current  density  occurs  under  the 
current  collector  land  areas  as  a  result  of  the  dominant  influ¬ 
ence  of  ohmic  losses  over  concentration  losses.  It  shows  that  the 
width  and  distribution  of  gas  channel  and  current  collector  land 
are  key  optimization  parameters  for  better  cell  performance.  The 
temperature  maxima  locate  in  the  cathode  CL  and  the  tempera¬ 
ture  variation  across  the  fuel  cell  increases  with  the  increasing 
of  current  density. 

Appendix  A 

According  to  the  experiment  of  real  gases,  the  thermal  con¬ 
ductivity,  viscosity  and  heat  capacity  of  each  gas  species  can  be 
described  by  polynomial  expression  of  temperature 

=  A0  +  AiT  +  A2T2  +  A3T\  (A.l) 

where  the  value  of  A/,  i  =  0,  . . .,  3  can  be  found  in  the  following 
Table  A.l. 
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